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LIKELIHOOD INFERENCE FOR PARTICLE LOCATION IN 
FLUORESCENCE MICROSCOPY 1 

By John Hughes, John Fricks and William Hancock 

Pennsylvania State University 

We introduce a procedure to automatically count and locate the 
fluorescent particles in a microscopy image. Our procedure employs 
an approximate likelihood estimator derived from a Poisson random 
field model for photon emission. Estimates of standard errors are 
generated for each image along with the parameter estimates, and the 
number of particles in the image is determined using an information 
criterion and likelihood ratio tests. Realistic simulations show that 
our procedure is robust and that it leads to accurate estimates, both 
of parameters and of standard errors. This approach improves on 
previous ad hoc least squares procedures by giving a more explicit 
stochastic model for certain fluorescence images and by employing a 
consistent framework for analysis. 

1. Introduction. The accurate and precise tracking of microscopic fluo- 
rescent particles attached to biological specimens (e.g., organelles, membrane 
proteins, molecular motors) can give insights into the nanoscale function and 
dynamics of those specimens. This tracking is accomplished by analyzing 
digital images produced by a CCD (charge-coupled device) camera attached 
to a microscope used to observe the specimens repeatedly. In this paper we 
introduce an improved technique for analyzing such images over time. Our 
method, which applies maximum likelihood principles, improves the fit to 
the data, derives accurate standard errors from the data with minimal com- 
putation, and uses model-selection criteria to "count" the fluorophores in 
an image. The ability to automate the process and quickly derive standard 
errors should allow for the analysis of thousands of images obtained from 
a typical experiment and aid in methods to track individual fluorophores 
across sequential images. 
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In fluorescence microscopy, a specimen of interest is tagged with a flu- 
orescent molecule or particle. The fluorescence microscope then irradiates 
the specimen with light at the excitation wavelength of the fluorophore, and 
when the excited electrons revert to the ground state they emit photons at 
the emission wavelength. A filter separates the emitted light from the exci- 
tation light so that only the light from the fluorescent material can pass to 
the microscope's eyepiece and camera system [Rost (1992)]. 

In general, the Rayleigh criterion implies that the maximum resolution 
for a light microscope should be roughly 250 nm (half of the wavelength 
of visible light); however, Selvin and his collaborators found that by fitting 
the center point of the point spread function one can locate a particle of 
interest. This technique is known as FIONA (Fluorescence Imaging with 
One-Nanometer Accuracy) and was introduced in Yildiz et al. (2003). The 
key element of FIONA is to focus attention on single fluorophores used as 
markers in biological specimens [Kural, Balci and Selvin (2005)]. By ana- 
lyzing sequences of images, molecular motors (e.g., myosin VI and kinesin) 
and other specimens can be tracked through time, giving researchers insight 
into their dynamics and biological function. For instance, Yildiz et al. used 
FIONA to find compelling support for the hypothesis that myosin V walks 
hand over hand and evidence to eliminate other hypotheses [Yildiz et al. 
(2003); Kural, Balci and Selvin (2005)]. 

A number of analysis techniques have been proposed for FIONA 
images. In 2001, Cheezum, Walker and Guilford compared four 
methods — cross-correlation, sum-absolute difference, centroid, and Gaus- 
sian fit — and ultimately recommended the Gaussian-fit method for single- 
fluorophore tracking. In the Gaussian-fit approach, the method of ordinary 
least squares (OLS) is used to fit a sum of symmetric bivariate Gaussian 
functions to the image. Least squares fitting is relatively efficient, and soft- 
ware to do it is widely available. Thompson, Larson and Webb subsequently 
proposed a "Gaussian mask" algorithm that is easier to implement than the 
Gaussian-fit method, is computationally less intensive, and performs nearly 
as well in simulations (2002). The Gaussian-mask algorithm is essentially a 
centroid calculation that weights each pixel with the number of photons in 
the pixel and with a bivariate Gaussian function integrated over the pixel. 
In both cases, simulation studies using typical experimental values showed 
that sub-pixel or even nanometer resolution was possible. 

The above mentioned Gaussian-fit and Gaussian-mask methods, while 
appealing, share two shortcomings. Since one or more beads may move out 
of frame for a particular image, the number of beads from one image to 
the next is not known a priori and must be determined for each image. 
Previous authors have attempted to solve this problem by means of a grid 
search, the first step of which is to scan the image for all pixels greater 
than some arbitrary threshold value. Each of these extreme pixels is taken 
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to be a bead location, and some region surrounding each extreme pixel is 
extracted from the image and processed by OLS or the Gaussian mask. 
Thompson et al. suggest a threshold that is eight standard deviations above 
the mean pixel value, but no explicit evidence is given in support of this 
choice [Thompson, Larson and Webb (2002)]. The correct threshold level 
for a set of images could possibly be approximated using simulations, but 
this would be a complex and computationally intensive task that would be 
necessary for each set of images, since the level of background noise may 
vary significantly from one experiment to the next. 

The second drawback is that estimates of precision are derived from simu- 
lation studies alone. If the probability model for the problem is misspecified, 
then error estimates based on simulations from that model will be inaccurate 
even if reasonable parameter values were used in the studies. And these val- 
ues will vary from image to image due to changing experimental conditions, 
for example, elevated background noise or slight changes in focus. A possi- 
ble solution is to perform a Monte Carlo simulation study using parameter 
values derived from the current experiment. But, given that the fitting pro- 
cedures themselves are time consuming, these approaches to standard error 
calculation may prove infeasible. It takes our algorithm several minutes — on 
a dual 2.8GHz Quad-Core Intel Xeon Mac Pro — to process an image with 
fifteen particles, which implies that bootstrapping standard errors for such 
an image would require hours of computation. Moreover, a full analysis of 
an experiment requires processing many hundreds of images. 

In what follows we present a new approach to counting and locating fluo- 
rophores. Our approach eliminates the need for a grid search and estimates 
standard errors from the data, without additional simulation, via standard 
likelihood tools. In Section 2 we present an explicit probability model for a 
FIONA image along with a maximum likelihood estimation procedure suit- 
able for this model. In Section 3 we discuss the properties of the approximate 
likelihood estimator presented in Section 2. In Section 4 we discuss stepwise 
model selection, which allows our procedure to automatically determine the 
number of beads in an image in a consistent manner. In Section 5 we describe 
the results of realistic simulation studies that support the approximations 
presented in Section 2 and demonstrate the robustness of our procedure. 
Finally, in Section 6 we carry out a complete analysis of an experimentally 
collected FIONA image, introducing relevant diagnostic criteria for our fit. 

2. Model for data from a single FIONA image. We first develop a model 
for the photon emission from particles distributed over a microscope slide. 
Photon emission from a constant source generally follows a Poisson distri- 
bution; this fact naturally leads to a Poisson random field model for the 
emission from a slide. We then express the effect of pixelation on the field, 
representing the emission as viewed from the digital camera, and employ 
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a normal approximation to the Poisson distribution. We arrive at our final 
approximate model by accounting for additional error introduced by the 
camera and its associated equipment. 

We begin with the standard model for the photon-emission pattern (as 
distorted by the point-spread function of a microscope objective) of a collec- 
tion of fluorophores distributed at random over some region of M 2 [Cheezum, 
Walker and Guilford (2001); Thompson, Larson and Webb (2002)]. 

Let N, a Poisson random field on a rectangular subset T of M 2 , represent 
the emission pattern of the sample. The intensity function can be defined 
for any Borel set R C T as 

(1) E{N(R)} = J I r \b + f> 3 • exp(- {X ~ X JI± {V ~ ^ ) } dxdy. 

Thus, N(R) is a Poisson random variable with the mean equal to a sum of J 
Gaussian functions, one for each bead, with Gaussian function j symmetric 
about (xj,yj) (which is contained in T). In addition, there is a constant 
background intensity of magnitude B representing background fluorescence. 
Although the intensity function for a bead is more often modeled in the 
physics literature by an Airy function, a Gaussian function approximates 
the Airy function quite well, and so we take the Gaussian centered at (xj,yj) 
to represent the (distorted) emission of bead j [Saxton (1997); Thompson, 
Larson and Webb (2002)]. 

The photons emitted by the sample are collected by a camera, the pixels 
of which can be represented by partitioning T into a uniform grid, where 
each pixel in the grid is square with side length (a) nm. Then, for a given 
pixel Z{ with center (xi,yi), 

E(Zi) 
(2) 



r-yi+a/2 


r-Xi+a/2 ( 


Jyi-a/2 . 


lxi-a/2 { 




( (x-xj) 2 + {y-yjf 



which is approximately equal to 

(3) (fl + Ev«g(- ( *~^ te ~ W) ' U«'- 



J 



Since a 2 is a constant, we allow the Aj and B to absorb it, arriving at 
(4) E(Z,) , B + E ^ ■ exp(- (I ' - - »> 2 ) = f. 
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Moreover, B is generally large enough to justify using a normal approxi- 
mation to the Poisson distribution. More precisely, 

(5) Zi~N(fi,fi), 

where we note that /j depends on the parameters of interest. In this model 
it is obviously important that B and the Aj's be constrained so that fi is 
nonnegative. 

The discretized Poisson random field described above is taken as the un- 
derlying model for the photon emission; however, additional error, which 
we will call instrumentation error, arises from various sources such as sig- 
nal quantization and dark current, an electric current that flows through 
a CCD even when no light is entering the device [Bobroff (1986); Thomp- 
son, Larson and Webb (2002)]. If we model the instrumentation error as a 
iV(0, 6) random variable independent of the intensity, then we have as a final 
approximate model for the data 

(6) Zi~N(fi,fi + 6). 

Consequently, the approximate likelihood of a given image with n pixels 

is 

m ^=n^j-p(-^). 

where /3 = (xq, yo, Ao, . . . , yj-i, -Aj-i, S, B, 9) T , the parameters of in- 
terest. This implies that the log-likelihood, without unnecessary constants, 
is 

(8) W = - $>(/, + 

i i 

which we maximize with respect to (3 to obtain /3mle- 



3. Estimating standard errors. From the theory of maximum likelihood 
estimators we know that, provided certain regularity conditions are met, 
a properly scaled MLE converges asymptotically to a normally distributed 
random variable. In our case, 

(9) [x n (p)] 1/2 &n-P) Ajv p (o,y, 

where X n (j3) is the Fisher information about (3 contained in a sample of size 
n, p is the dimension of (3, and I is the p x p identity matrix. This implies 
that 
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for large n, and so the diagonal elements of [X n (/3)] 1 are approximate sam- 
pling variances for the estimators (3 n . However, we do not know [X n (/3)]~ 1 
because the true f3 is unknown and analytical calculation of the information 
is prohibitively complicated. 

Consequently, we use the standard substitution [X n (/3)] _1 using the ob- 
served information, J n (P n ) = i~ df^dp/ niPn)], tliat is > [ X n{P)]~ l ~ 

[J'n{P n )]~ 1 ■ Estimating [X n (/3)] _1 by inverting t 7 n ( / 9 n ) has the advantage 
of avoiding closed form derivatives which are unwieldy in this case. 

The preceding implies that the joint distribution of Xj and y 3 - is approxi- 
mately bivariate normal. More precisely, 



(11) £j = (.Xj,yj) T ~ N 2 [Hj = {x j ,y j ) T ,Y, 



oh 

where 

(12) £ ■ = [JniPnT 1 = \ Covixj^j) 

_Cov{xj,yj) Vax(yj) 

Since the contours, that is, the equidensity curves, of the bivariate nor- 
mal distribution are ellipses, an approximate 95% confidence region for the 
location of bead j also takes the form of an ellipse: 

(13) {»j ■ (Aj - Vjf^j ~ = Xo.95,2}> 

where Xo 95 2 denotes the 95th percentile of the x 2 distribution with 2 degrees 
of freedom [Ravishanker and Dey (2002)]. For an image with multiple beads, 
the typical case, we may want a collection of random ellipses, the union of 
which will enclose all J beads with probability 0.95. We can accomplish 
this by using the Bonferroni correction, which assigns to each bead an error 
rate of thereby making the image wide error rate 0.05. The resulting 
collection of simultaneous confidence ellipses is given by 

(1 4 ) {Vj ■ (Aj - Vjf^j (Aj - = X?-0.05/J,2}- 

We evaluated our standard-error estimation and the convergence of our 
estimator by way of a simulation study. Ten thousand 100 x 100 single- 
bead images were simulated, each image having its lone bead located at 
(7823,3353), where the coordinates are given in nanometers from the lower 
left corner. (For an idea about the nature of the data, see Figure 2, which 
has four beads.) Table 1 shows /3mle f° r a single image along with standard- 
error estimates for that image and the true standard errors gleaned from all 
10,000 images. Our estimated standard errors are in close agreement with 
the true standard errors. 

Figure 1 shows estimated densities for the sampling distributions of xq, 
i/0) Aq, S, B, and 9, respectively. Appropriate normal densities (dashed) are 
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shown superimposed. Each normal density is centered at the true value for its 
parameter. It is clear that the sampling distributions converge to normality, 
but the estimators of S and Aq are slightly biased in opposing directions; 
intuitively, estimation of S works in opposition to that of Aq. This is because 
the fitting procedure is attempting to simultaneously conform S to the base 
of the Gaussian peak and Aq to the peak's height, and an adjustment of 
either estimate nudges the other in the opposite direction. 

4. Model selection: how many beads are in the image. Previous authors 
have suggested that the number of beads in an image be determined by ap- 
plying a grid search algorithm prior to fitting [Cheezum, Walker and Guilford 
(2001); Thompson, Larson and Webb (2002)]. As we mentioned in Section 1, 
any pixel with an intensity above some threshold is identified with a bead, 
and then some region in the vicinity of the pixel is extracted from the im- 
age and fitted. This thresholding approach may be adequate for producing 
initial estimates of bead locations, but thresholding prevents full automa- 
tion because the threshold must be chosen by the investigator. And even a 
seemingly well-chosen threshold may be too large to distinguish dim beads 
from background noise. 

Our procedure eliminates these problems by fitting first and selecting 
the number of beads based on those fits. In our scheme, an approximate 
information criterion derived from an OLS fit is used to approximate the size 
of the model, and then, because the candidate models are nested, likelihood 
ratio tests are used to select the final model. As we will show in Section 
5, this approach is able to identify all of the beads, even very dim ones, 
automatically. 

Our algorithm has a preliminary stage for estimating the number of beads 
and producing initial estimates of all parameters except 9, and a final stage 
for estimating 9, giving maximum likelihood estimates of the other param- 
eters, and accurately determining the number of beads. The preliminary 



Table 1 

Estimation of standard errors. Estimation error is the 
true value of the parameter minus the estimated value 



Parameter 


Truth 


Est error 


Sim SE 


SE 


•f'o 


7823 


-0.4 


0.420 


0.422 


yo 


3353 


0.1 


0.424 


0.421 




15,000 


99.8 


77.4 


61.2 


s 


200 


-0.2 


0.393 


0.341 


B 


200 


0.1 


0.173 


0.171 





100 


6.7 


4.26 


4.16 



<s 



J. HUGHES, J. FRICKS AND W. HANCOCK 





Fig. 1. These plots show density estimates for the sampling distributions of xo, yo, Aq, 
S, B, and 8, respectively, with normal densities (dashed) superimposed. 
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stage assumes zero beads initially and fits f(x,y) = B to the image using 
OLS. Using the least squares fit at each stage, the information criterion 

(15) IC ^= n \n(?^-\+ P ^i 

is computed, where k is the (assumed) number of beads, n is the sample size, 
RSS is the residual sum of squares, and p is the number of free parameters. 
Note that IC is an increasing function of RSS and p, which implies that IC 
rewards a better fit (smaller RSS) and penalizes more free parameters. On 
the next iteration, one bead is assumed, and so 

(x-x ) 2 + (y-y ) 2 ' 



(16) f(x,y) = B + A -exp 



S 2 



is fit to the image, producing IC^ l \ Iteration continues until IC^ k ' > IC^" 1 ^, 
which indicates that the image contains k — 1 beads. 

Note that IC is a nonstandard information criterion. We found that even 
the Bayesian Information Criteria (BIC) does not penalize additional pa- 
rameters sufficiently. This can allow the initial stage of the algorithm to 
significantly overestimate the correct number of beads, causing much unnec- 
essary computation during the final stage of the algorithm. Our simulations 
showed that replacing BIC's ln(ra) with y/n minimizes overfitting. 

As the algorithm makes an initial forward sweep over possible models, 
the OLS parameter estimates are saved. This initial sweep stops based on 
the information criteria, IC. Those estimates are used to initialize the max- 
imum likelihood estimation carried out in the final backward sweep which 
terminates using likelihood ratio criteria. Providing the MLE routine with 
good initial estimates of all parameters except 9 allows the MLE to converge 
faster than it otherwise would. 

The parameter fits and the selection criteria initially computed are then 
used to find the final parameter estimates and make the final model selection. 
The key differences are that OLS is replaced by MLE and IC is replaced by 
the likelihood ratio statistic 

(17) G\ beads) = -2{i n 0^) - i n 0^ +1) )h 

which should be approximately y 2 distributed with three degrees of freedom 
(because each bead is associated with three parameters: Aj, Xj, and yj). The 
full algorithm is given in pseudocode below. 

5. Simulated examples. In this section we present a series of simulated 
examples. We first apply our procedure to a typical image simulated from 
the Poisson plus Gaussian model presented above. Then we examine the 
robustness of our procedure by applying it in three atypical scenarios: dim 
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Algorithm 4.1 LocateBeads (pixels) 



beads 

0OLS S \ Icibeads) } <~ OLS(pixels, beads) 
repeat 

beads <— beads + 1 

0o£t\ JC( Ws )} <- OLS(pixe/s, beads) 
until ic^ beads ^ > ic^ beads ^ 
beads <— beads — 1 

'•(beads) "(beads) "(beads) 

Ipmle .4(Pmle )} ^MLE (pzxe/s, beads, /3 ols ) 
repeat 

beads <— 6eads — 1 

r "(beads) . « (beads). -. , "(beads). 

{Pmle -4(Pmle )}<-MLE(^xe<s, beads, /3 OLS J 

Gfbeads)* 2{£ n (/3 MLE ) — ^n(/^MLE )} 

UntilG ( & ea^)> 7 - 81 = xi 9 5,3 
"(f>eac2s+l) 

return (/3 M LE ) 



beads, beads in close proximity, and beads that are not entirely contained 
by the image. Finally, we investigate the sensitivity of our procedure to 
misspecification of the instrumentation error. 

First, we fit an image with fifteen roughly even-spaced fluorophores to ver- 
ify that our method can handle the substantial numbers that are sometimes 
found in experimental data. The parameter estimates and their approximate 
standard errors appear in Table 2. While the time to fit such a larger exam- 
ple is considerable, the method works well and finds the correct number of 
beads without difficulties. 

Figure 2 shows an image that contains a bead that is very dim (A = 
400) relative to the image's other three beads (A = 15,000). We simulated 
1000 such images and applied our procedure to each. Our procedure was 
able to estimate the bead's location to within a standard error of less than 
six nm, which, relative to the other beads, represents a tenfold decrease in 
resolution for a fortyfold decrease in brightness. Table 3 gives the results 
for this simulation study. Additionally, simulations showed our algorithm 
capable of consistently locating (against a background of 200) beads as dim 
as A = 75, which implies a contrast ratio, that is, the ratio of the brightest 
pixel value and the background value, equal to 1.4 (versus 75 for a typical 
bead). 

Figure 3 shows an image with two beads whose centers are separated by 
only 400 nm. We again applied our procedure to 1000 like images, each image 
having four beads. Our algorithm was able to distinguish the two close beads 
with only a slight loss of precision in the direction of the line between the 
beads, as is shown in Table 4. 
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Table 2 

Localization for a Typical FIONA Image. 
A = 15,000 for each bead, and S = 200, B = 200, 
6 = 100. No estimation error for A was greater 
than 105, and every approximate confidence 
interval save one covered the truth. The 
estimation errors for S, B, and 8 were 0.1, —0.1, 
and 1.7, respectively, and their confidence 
intervals covered the true values 



Parameter 


Truth 


Est error 


SE 


Xq 


23,566 


0.9 


0.422 


V/n 


4852 


0.9 


0.423 


X\ 


2522 


-0.6 


0.423 


//i 

U 1 


18,672 


0.1 


0.421 




10,475 


-0.6 


0.424 




4858 


—0.05 


0.423 


X3 


16,643 


-0.3 


0.422 


2/3 


19,505 


0.6 


0.423 


X4 


6842 


0.4 


0.423 




16,060 


0.1 


0.421 


■IT, 


17,753 


-0.3 


0.421 


ys 


28,518 


1 


0.423 


xe 


28,956 


0.2 


0.421 


,'/(. 


6771 


0.2 


0.421 


X7 


27,512 


0.3 


0.419 


2/7 


3454 


-0.4 


0.419 


Xg 


4165 


0.3 


0.422 


2/8 


13,466 


-0.5 


0.422 


Xg 


28,960 


0.4 


0.421 


V9 


11,712 


-0.5 


0.421 


Xio 


25,394 


0.1 


0.422 


llw 


28,468 


-0.4 


0.421 


X-ix 


29,112 


0.1 


0.423 


!)n 


28,770 


-0.6 


0.422 


Xl2 


18,028 





0.422 


2/12 


21,796 


0.1 


0.422 


Xl3 


28,318 


0.2 


0.434 


V13 


12,251 


0.2 


0.428 


X14, 


27,757 


0.8 


0.432 


2/14 


11,937 


0.5 


0.426 



The second image in Figure 3 shows a bead whose center is only 50 nm 
from the image's edge. Our algorithm was able to localize such beads with 
a loss of precision in the y direction that is quite acceptable and perhaps 
even surprisingly small given that nearly half of the bead is missing. Table 5 
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Fig. 2. An image with four beads, one very dim. Our procedure locates all four beads. 

reports the results for a 1000-image study, where again each image had four 
beads. 

Maximum likelihood estimation is often sensitive to model misspecifi- 
cation, and so we investigate the performance of our procedure when the 
instrumentation error is not N(0,6). The instrumentation error for each im- 
age has mean zero and variance 8, but otherwise the errors are distributed 
rather differently. The model was simulated, but with heavy-tailed (£3 dis- 
tributed) and asymmetric (exponentially distributed) instrumentation error. 
More specifically, the model was simulated according to 

(18) Zt ~ Poi(/0 + y/e/3t 3 
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Table 3 
Localization of a dim bead 



Parameter 


Truth 


Est error 


Sim SE 


SE 


Xq 


4021 


0.03 


0.431 


0.422 


yo 


5172 


-0.13 


0.440 


0.422 


id 


i ^ nnn 


RQ Q 
uy.y 


uo.ooy 




Xl 


1497 


-0.49 


0.430 


0.422 


m 


9241 


-0.04 


0.428 


0.422 


A, 


15,000 


-8 


61.846 


49.508 


X2 


7920 


-0.61 


0.425 


0.423 


3/2 


1807 


0.29 


0.405 


0.424 


A 2 


15,000 


—7.3 


60.757 


49.513 


x 3 


6000 


4.69 


5.10 


5.03 




8722 


11.89 


5.30 


5.08 


A :i 


400 


—8.761 


12.377 


11.11 


S 


200 


-0.174 


0.220 


0.197 


B 


200 


-0.194 


0.176 


0.174 





100 


2.575 


4.076 


4.25 






Table 4 






Loca 


lization c 


if beads in close 


proximity 




Parameter 


Truth 


Est error 


Sim SE 


SE 


Xo 


4021 


-0.22 


0.430 


0.422 


yo 


5172 


0.45 


0.418 


0.421 


A 


15,000 


22.5 


66.827 


54.260 


■n 


1497 


0.06 


0.424 


0.422 


yi 


9241 


0.36 


0.417 


0.423 


Ax 


15,000 


30 


62.375 


48.019 


X2 


7920 


-0.21 


0.467 


0.445 


V2 


1807 


1.89 


0.579 


0.555 


A 2 


15,000 


-17.5 


61.06 


48.328 


X3 


7920 


0.06 


0.449 


0.445 


in 


2207 


-0.85 


0.563 


0.552 


A, 


15,000 


13.8 


67.593 


54.291 


S 


200 


-0.119 


0.190 


0.178 


B 


200 


-0.196 


0.169 


0.173 





100 


6.518 


4.261 


4.182 



and according to 

(19) Zi ~ Poi(/i) + E W (Ve) - Ve, 

where Poi(A) denotes the Poisson distribution with rate A, t v denotes the t 
distribution with v degrees of freedom, and Exp (A) denotes the exponential 
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D O 0.2 D4 Q 6 0.8 1.0 




0.0 02 01 OB 1.0 

Fig. 3. The upper image shows two beads in close proximity. The lower is an image with 
a partial bead. 

distribution with mean A. Tables 6 and 7 show that localization was not 
affected by these misspecifications. Again, 1000 images were used for each 
study. 

Table 8 gives the coverage rates of our approximate 95% confidence regions 
for all of the previously mentioned simulation studies and for (first row) a 
study of 1000 typical images. The coverage rates clearly suffer a bit for all 
except the typical and asymmetric scenarios. 
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6. Analysis of an experimentally observed FIONA image. In this section 
we apply our procedure to an experimentally observed FIONA image shown 
in Figure 4. Table 9 shows our parameter estimates for this image. 



Table 5 
Localization of a partial bead 



Parameter 


Truth 


Est error 


Sim SE 


SE 


•en 


4021 


-0.09 


0.413 


0.424 


2/0 


5172 


-0.03 


0.410 


0.423 


A 


15,000 


143 


69.09 


48.217 


•Cl 


1497 


-0.39 


0.421 


0.424 


Vi 


9241 


0.39 


0.432 


0.424 


A, 


15,000 


143 


69.581 


71.17 


■f-2 


7920 


0.37 


0.430 


0.422 


2/2 


1807 


0.35 


0.417 


0.422 


A 2 


15,000 


104.4 


64.671 


48.416 




6000 


0.11 


0.524 


0.525 


2/3 


50 


-0.822 


0.924 


0.846 


A, 


15,000 


37.2 


62.715 


48.738 


S 


200 


-0.51 


0.214 


0.187 


B 


200 


-0.119 


0.169 


0.175 





100 


0.588 


4.243 


4.265 



Table 6 

Estimation when instrumentation error is heavy-tailed 



Parameter 


Truth 


Est error 


Sim SE 


SE 


Xq 


4021 


-0.25 


0.433 


0.422 


2/o 


5172 


-0.29 


0.409 


0.422 


An 


15,000 


-18.3 


59.564 


47.923 


XI 


1497 


-0.3 


0.437 


0.421 


Hi 


9241 


-0.54 


0.425 


0.421 


Ai 


15,000 


26.9 


65.228 


47.658 


x% 


7920 


-0.19 


0.428 


0.422 


2/2 


1807 


-0.13 


0.422 


0.422 


A 2 


15,000 


33 


58.357 


47.671 


X3 


6000 


-0.6 


0.432 


0.421 


2/3 


8722 


-0.1 


0.413 


0.422 


A, 


15,000 


-19.8 


65.419 


47.830 


S 


200 


0.073 


0.190 


0.171 


B 


200 


-0.017 


0.175 


0.173 





100 


6.435 


9.759 


4.191 
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To evaluate the fit to the image, we ran numerous diagnostics to verify 
that the observed data originates from our random field model. Our approx- 
imate model implies that for the ith pixel 

(20) Zi~N(fiJi + 0), 

spatially independent of the other pixels. If our model is correct, then we 
should have, for the ith. error, 

(21) Ei = ^A~iV(0,l), 

spatially independent of the other errors. This implies that the variogram 

(22) 7 (h) = ±Var{e(s + h)-e(s)} 

should equal one for all locations s and lag (displacement) vectors h. 

We plot empirical (residual) variograms to determine spatial indepen- 
dence and use a normal probability plot to check for normality [Cressie and 
Hawkins (1980)]. The plots for our example image are shown in Figure 5. 
Except for some anomalous features in the lower tail of the probability plots, 
the diagnostics give a good indication that our proposed model is sound. The 
standardized residuals were also checked and no blatant violations of what 
would be expected for independent, identically distributed data were found. 



Table 7 

Estimation when instrumentation error is asymmetric 



Parameter 


Truth 


Est error 


Sim SE 


SE 


Xo 


4021 


-0.26 


0.428 


0.420 


yo 


5172 


-0.75 


0.398 


0.421 


A 


15,000 


5.1 


60.595 


47.619 


XI 


1497 


-0.39 


0.410 


0.421 


J/1 


9241 


0.07 


0.422 


0.422 


Ai 


15,000 


-149.7 


62.130 


48.146 


X2 


7920 


-0.17 


0.413 


0.421 


2/2 


1807 


0.08 


0.442 


0.421 


A 2 


15,000 


-68 


57.774 


47.940 


X3 


6000 


0.67 


0.424 


0.420 


V3 


8722 


0.04 


0.420 


0.419 


A, 


15,000 


-74 


61.813 


48.109 


S 


200 


0.12 


0.188 


0.170 


B 


200 


0.44 


0.178 


0.173 





100 


5.075 


4.798 


4.203 
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Table 8 

Coverage rates of approximate 95% confidence ellipses 



Scenario 


Bead type 


Coverage rate 


Typical 


typical 


n A AG/ 

94.470 




typical 


95.6% 




typical 


95.1% 




typical 


95.4% 


Uim 


typical 


yo.D/o 




typical 


94.6% 




typical 


95.2% 




dim 


93.7% 


Ulose 


typical 


9O.1/0 




typical 


94.6% 




close 


94.3% 




close 


93.8% 


Partial 


typical 


95.3% 




typical 


94.5% 




typical 


95.8% 




partial 


93.0% 


Heavy-tailed 


typical 


94.9% 




typical 


94.3% 




typical 


94.8% 




typical 


93.5% 


Asymmetric 


typical 


96.1% 




typical 


95.6% 




typical 


95.8% 




typical 


95.2% 




Table 9 




Parameter estimates for 


an 


experimentally observed FIONA image 


Parameter 


Estimate 


SE 


Xq 


12,168.1 


4.83 


yo 


7570.4 


4.79 


A 


269.8 


10.5 


X\ 


12,296.1 


4.71 


yi 


16,509 


4.63 


A! 


275.7 


10.3 


S 


175.7 


3.31 


B 


33.4 


0.0418 





80.8 


0.632 
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0.!t Z U.4 0.6 0.8 1.0 




Fig. 4. An experimentally observed FIONA image. 

7. Conclusions. The method outlined in this paper allows for the auto- 
mated analysis of FIONA images, including the ability to select the number 
of fluorophores in an image. By using a likelihood framework, the method 
also allows for standard errors to be calculated simultaneously with the es- 
timates. The method was then verified through simulation and the analysis 
of collected data. We hope that this case study will serve as an example of 
applying traditional statistical theory to enhance the analysis of nanoscale 
experimental methods where algorithmic approaches have been favored. 

Since this method is largely automated through model selection tech- 
niques, it can handle the analysis of "movies" by processing each frame. 
Since the method also returns standard errors for the locations of the flu- 
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Fig. 5. A nanogram plot and a normal probability plot of the standardized residuals for 
the real FIONA image. The lags for the variogram plot are given in nanometers, and the 
image in question was approximately 27,000 nm on a side. 

orophores, this opens the possibility of creating tracking methods to follow 
dynamic specimens using not only the position data but the information on 
observational errors which are given. 

Estimation and model selection can be done using our free C++ applica- 
tion, beads, which makes extensive use of the GNU Scientific Library, and fit 
diagnostics can be carried out using our free R software package, FIONAdiag 
[Galassi et al. (2008)]. 
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